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We introduce a patch to the commonly used public codes CAMB and CosmoMC that allows the 
user to implement a general modification of the equations describing the growth of cosmological 
perturbations, while preserving the covariant conservation of the energy-momentum. This patch 
replaces the previously publicly released code MGCAMB, while also extending it in several ways. 
The new version removes the limitation of late-time-only modifications to the perturbed Einstein 
equations, and includes several parametrization introduced in the literature. To demonstrate the 
use of the patch, we obtain joint constraints on the neutrino mass and parameters of a scalar-tensor 
gravity model from CMB, SNe and ISW data as measured from the correlation of CMB with large 
scale structure. 

I. INTRODUCTION 

Future cosmological surveys will have the ability to measure the growth of large-scale structure with accuracy 
sufficient for discriminating between different models of dark energy (DE) and modified gravity (MG). To this aim, 
an interesting approach is to parametrize a general theory of gravity in a model independent way so that it reduces 
to general relativity (GR) for a choice of parameter values, and can also be related to certain classes of MG models. 
This approach has been taken by several groups to test the standard Lambda Cold Dark Matter (ACDM) model and 
to constrain parameters of various DE and MG models [1-16]. We refer the reader to [17, 18] for a comprehensive 
review of the recent advances in the field of MG and ways of testing GR on cosmological scales. 

Code for Anisotropics in the Microwave Background (CAMB) [19, 20] uses the equations of GR together with 
Boltzmann (conservation) equations to calculate cosmological observables. It needs to be modified for studying the 
alternative gravity theories, as for example was done in [3, 5, 6, 12]. One then needs a system of equations that is 
meaningful across a wide range of scales and redshifts. As stressed in [3], the evolution of perturbations on super- 
horizon scales is governed by a set of consistency conditions which are separate from the sub-horizon dynamics. Namely, 
as shown in [21, 22], in the absence of entropy perturbations, the space curvature defined on hyper-surfaces of uniform 
matter density, £, must be conserved on scales outside the horizon in order to be consistent with the overall expansion 
of the universe. Hence, a consistent system of equations should decouple the super- and sub-horizon regimes. This 
separation of scales is made explicit in the parametrized Post-Friedmannian (PPF) Framework of [3], where a different 
systems of equations are used on super-horizon and sub-horizon scales. The advantage of the method advocated in 
this paper and used in MGCAMB (Modification of Growth with CAMB) [6] is that it employs a single system of 
equations across all linear scales, without sacrificing any of the important consistency conditions. The super-horizon 
and sub-horizon evolution decouple naturally, without having to be explicitly separated. 

Here we introduce a generalized version of MGCAMB, which now has a form of a patch to CAMB and CosmoMC, 
available for download and public use at [23] . This version is capable of evolving all the species in the universe from 
early times to present. While the previous version was based on CAMBsources [24], which made it difficult to use in 
CosmoMC, this version is based on CAMB and can be used easily with CosmoMC to constrain desired parameters. 
Several parametrizations from the literature are included in this version which makes it useful for investigating a wide 
range of models. 

Massive neutrinos are not part of the standard ACDM model, and can be considered as an extension of the model 
that modifies the growth of large scale structure. It is, therefore, interesting to study possible degeneracies between 
effects of the neutrino mass and various MG or dark energy models which also modify the growth on large scales. We 
used MGCAMB to study the degeneracy between massive neutrinos and a specific class of f(R) and Chameleon- type 
models using a combination of CMB, supernovae and the CMB-LSS cross-correlation data of Ho et al. [25]. We find 
the degeneracy to be small due to the smallness of the overlap of the scales at which the respective modifications 
affect the observables we have considered. 

This paper is organized as follows. In section II, we review the linear scalar perturbation equations for GR and 
introduce our modification to these equations. Some of the alternative parametrizations in the literature are also 
reviewed in this section. In Section III, we describe the implementation of our equations in CAMB and introduce 
the new MGCAMB patch. Joint constraints on massive neutrinos and some MG models are obtained in Section 
IV. Certain details of the way the modification was implemented in the CMB source function in CAMB and in the 
CMB-LSS cross-correlation patch of Ho et al. [25] are given in two appendices. We summarize in Section V. 
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II. FORMULATION 



As a starting point we assume that the background universe is described by a flat Friedmann-Robertson- Walker 
(FRW) metric, and consider linear perturbations in the metric and the energy-momentum tensor. In the Newtonian 
gauge, the perturbed line element has the form 

ds 2 = a(r) 2 [-(l + 2^)dr 2 + (1 - 2<£)dx 2 } , (1) 

where r is the conformal time and vff and $ are the two scalar metric potentials. In the synchronous gauge, which is 
used in CAMB, the perturbed line element is 

ds 2 = a{T) 2 [-dr 2 + (Sij + hi^dx'dx*] , (2) 

where the scalar part of hij is represented by functions rj and h which in Fourier space are defined as [26] 

h i:j (x, t) = J d 3 ke l ^ s ^kikjh(k, r) + (kkj - 67 ? (fc, r)| , k = kk . (3) 

We also consider perturbations in the radiation-matter content, which includes photons, neutrinos, baryons and cold 
dark matter (CDM): 

n + ST° = -p(l + 6) , 
Tf + 6T? = -(p + P)vi, 

T] + <5T/ = (P + 5P)S) + ttJ , (4) 

where S = Sp/p is the density contrast, v is the velocity field, 5P is the pressure perturbation and 7rj denotes the 
tracclcss component of the energy-momentum tensor perturbations. Unless noted otherwise, all the quantities are 
sums over all the species present in the universe. Working in Fourier space, it is convenient to introduce the comoving 
density perturbation A: 

P A = p6 + 3j(p + P)v , (5) 

where % = a/ a, and throughout this paper the overdot denotes the derivative with respect to the conformal time r. 
Note that we also use the physical Hubble parameter, H = H/a, later in the text. The anisotropic stress a and the 
momentum perturbation 6 arc defined as 

{p + P)a = -{k%- l -5^i , (6) 
(p + P)6 = ik j 6T9 . (7) 

The energy-momentum conservation (T.^ = 0) provides two equations relating the metric potentials and pertur- 
bations in each fluid. In the Newtonian gauge they are [26]: 

6 = -^(1 - 3w)6 - -^6 + k 2 6 - k 2 a + kH . (8) 

a 1+w 1 +w w 

where w = P/p. The anisotropic stress a vanishes for baryons and CDM, while for relativistic species, i.e. photons 
and neutrinos, it is generated through free-streaming and is related to 8 and 9 via Boltzmann equations. 

Two more equations are needed to close the system of equations for metric and energy- momentum equations. For 
example, in GR, the Poisson equation relates the comoving density perturbations to the metric potential <!>, while the 
anisotropy equation relates the two metric potentials and the anisotropic stress: 

k 2 <P = -47rGa 2 pA , (9) 
k 2 ($ - \H) = 12nGa 2 (p + P)a . (10) 
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In an alternative theory of gravity these relations would generally be different. To parametrize possible departures 
from the ACDM growth at late times {e.g. redshifts z < 30 when the contribution of relativistic species can be 
neglected) the following parametrization was used in the original version of MGCAMB [6]: 

k 2 q> = -4T:Ga 2 ti(k, a)pA , (11) 

i= 7 (M), (12) 

where the two scale- and time-dependent functions p(k,a) and 7(fc, a) were introduced to encode any modification 
to (9) and (10). Note that, in this formulation, it is the Newtonian potential ^ that appears on the right hand side 
of (11), while the curvature potential $ appears in the original Poisson equation (9). This is motivated by the fact 
that none of the observables depends solely on For example, weak lensing probes the combination ($ + "J), while 
clustering of matter and the peculiar velocities are most directly related to $ , which follows from the conservation 
equations (8). Thus, measurements of peculiar velocities and galaxy counts (up to bias factors) can constrain the 
function p most directly. Still, there can be circumstances in which it is advantageous to use other related parametric 
forms, e.g. when making a connection with particular theories [30]. We discuss some of the alternatives in II B. 

It should be emphasized that /i and 7 do not necessarily have a simple form in specific models of MG, and generally 
depend on the choice of the initial conditions 1 [27-29]. For instance, in scalar-tensor models of gravity, the ratio of 
$ and ^ is not a fixed function of k and a. Instead, it is an expression that involves the time derivatives of Vl 7 and 
$. This means that p and 7 correspond to solutions of equations of motion of a theory, rather than being a general 
prediction of a theory. 

It is also important to make the distinction between the formalism, based on writing modified linearized Einstein 
equations in terms of functions fi and 7, and any particular parametric form that one may chose for these functions. The 
formalism is self-consistent, and satisfies the super-horizon consistency condition [21, 22] as long as (k / (aH)) 2 / (fij) —> 
in the k/{aH) — > limit [30]. It allows for model-independent tests of GR by simply looking for deviations of the 
functions from unity. If one wanted to use this formalism to compute observables in a specific theory one, strictly 
speaking, would have to first solve the equations of motion for perturbations in that theory, find /1 and 7 by taking 
the appropriate ratios of the solutions, then substitute them into MGCAMB. Another approach, which would involve 
a degree of approximation, would be to design parametric forms for fi and 7 which can mimic solutions of particular 
theories [4, 7, 31, 32]. One needs to be careful when interpreting the constraints on /1 and 7 in terms of theoretical 
parameters, paying attention to the choice of the initial conditions and other possible assumptions and simplifications. 
We discuss a couple of specific parametric forms in II A. 

Eqs. (11) and (12) are only applicable to late times (low redshifts) since the anisotropic stress term, which is 
important at epochs when radiation and relativistic neutrinos are relevant, is ignored. Here, we generalize these 
equations to be valid at all times: 

fc 2 * = -n(k, a)ATrGa 2 {pA + 3(p + P)a} , (13) 
fc 2 [$- 7 (fc,a)*] = f i(k,a)l2nGa 2 (p + P)a , (14) 

where, in Eq. (13), the anisotropic a is added so that Eqs. (9) and (10) are correctly recovered in the ^ = 7 = 1 limit. 
In this form, the equations include all the species in the universe and can be used for all times, assuming one uses 
appropriate initial conditions. In the current version of the MGCAMB, we kept the CAMB default initial conditions 
corresponding to 7 = 1 when a given k is well outside the horizon during the radiation era (p becomes irrelevant on 
superhorizon scales [30]). We plan to generalize the initial conditions to be consistent with a general choices of 7 in 
future upgrades of MGCAMB. Evolving Eqs. (13) and (14) with given functions p and 7, along with the conservation 
and Boltzmann equations, allows us to compute various cosmological observables that can then be compared to the 
ACDM results. 



A. Specific parametric forms of p and 7 

While MGCAMB works for any form of p(k, a) and j(k, a) supplied by the user, we have coded in several specific 
paramctrizations previously used in the literature. For example, a particular form of p and 7 that holds approximately 
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in some classes of f(R) and scalar-tensor theories was introduced in [4], hereafter referred to as the BZ parametrization: 

1 + j8iA? k 2 a s 



H{k,a) 



1 + Af k 2 a 
l + j3 2 X% k 2 a s 



7(fc ' Q)= 1 + AgPa* ' (15) 

where /3j's are dimensionlcss couplings and Aj's have dimension of length. The previous version of MGCAMB [6] has 
been used to calculate cosmological observables and to forecast errors on the five BZ parameters for future weak 
lensing surveys, such as DES and LSST, in combination with CMB data from Planck. 

Specific sets of BZ parameter values can be chosen to correspond to particular f(R) [31, 33-37] and Chameleon 
type [38, 39] models. Moreover, not all of the parameters are independent. For example, f(R) models with action 
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16ttG 



d 4 xV^g[R + f(R) + C m ] , (16) 



can be tuned to reproduce any background expansion history, and the remaining relevant quantity is the squared 
Compton wavelength of the new scalar degree of freedom Jr = df / dR mediating the fifth force. In units of the Hubble 
length squared it is given by [34, 37] 

Bs J**_ jmUj^H\-\ (17) 



1 + /r d In a \ d In 

where R is the background Ricci scalar. Thus, for a fixed background expansion history, different f(R) models can 
be parametrized by the parameter Bq, which is the value of B today. It was suggested in [7] that for Bq < 1 the large 
scale growth in f(R) models can be modelled using a BZ form with (3\ = 4/3, A| = PxX 2 and /3% — 1/2 as: 

1 1 + |Aj fc 2 a 4 

M(fc ' fl) ~ 1- 1.4- 10-8|Ai| 2 a3 1 + A 2 fc 2 a 4 > ( 18 ) 

where the pre-factor is introduced to account for the background rescaling of the Newton's constant which can be 
important in some models, A 2 = Bq c 2 /(2Hq), and 7 is given by the BZ form (15). We will use this parametrization 
in the next section to demonstrate the use of the updated MGCAMB and to constrain the Bo parameter. 

Similarly, models with a Yukawa-type dark matter interaction, e.g the Chameleon type models on large scales, 
can also be approximately expressed in terms of the extended BZ parametrization of [7]. In this type of models, the 
allowed ranges of the BZ parameter values arc < B = 2X\Hq/c 2 < 1, < /3i < 2 and 1 < s < 4, while X 2 = PiXf 
and fa — 2//?i — 1, where the parameter Bq is again related to the Compton wavelength of the extra scalar degree of 
freedom. As shown in [7], and in the next section, the data considered in this paper cannot constrain parameters Bq 
and s in this model. 

Parametric forms such as (15) or (18) are certainly not guaranteed to be accurate in representing solutions of all 
f(R), Chameleon or other MG models at all values of the theoretical parameters. An extensive investigation of the 
range of their applicability would be a worthwhile pursuit which, however, is outside the scope of this paper. The 
validity of the original BZ form (15) was questioned in [32], where a more accurate parametric form specific to f(R) 
was suggested. In Section IV we demonstrate the use of MGCAMB by working with the modified BZ form (18) of 
[7] that has the additional pre-factor in /1. This pre-factor becomes important at Bq ~ 1, which we confirmed by 
reproducing the CMB spectrum plots in [37, 42] for Bo < 2 where exact f(R) equations of motion were used 2 . We 
also note that our bounds on the f(R) parameter B are in agreement with [15], who used the PPF parametrization 
of f(R) solutions described in [40]. 

B. Alternative parametrizations 

The parametrized modification of perturbed Einstein equations defined by Eqs. (13) and (14) is not unique. Other 
forms may be useful when working with specific cosmological observables or when testing particular theories. Generally, 
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it should be possible to express these other parametrizations in terms of p and 7. For example, in certain circumstances, 
it can be desirable to work with functions p and S, defined as 

k 2 ^ = -4nGa 2 p(k, a)pA , 

k 2 (<P + ty) = %nGa 2 Y>(k 1 a)pA. (19) 

In [14] it was shown that working with these functions clearly captures the complementary information in weak lensing 
and peculiar velocity experiments. In the limit of negligible anisotropic stress er, the function £ is simply related to p 
and 7. This, and other similar possibilities were discussed in [30] and it should be easy to modify MGCAMB to work 
in such cases. 

In [12], functions Q(k,a) and R(k,a) were introduced as 

fc 2 $ = -4nGa 2 QpA, 

fc 2 (* - = -l2nGa 2 Q(p + P)a , (20) 

which we include as an alternative parametrization in MGCAMB. Namely, when the appropriate option is selected, 
Eqs. (20) are evolved in place of (13) and (14). Functions Q and R are simply related to p and 7 in the limit of 
negligible anisotropic stress a: 

Q = M7 , R = r X , (21) 

although we do not rely on this conversion in practice. As a consistency check, we have reproduced some of results in 
[12] using MGCAMB. 

Another commonly used parameter is the growth index 7, proposed in Ref. [41]. In this paper we will denote the 
growth index 7 as 7^ to distinguish it from our function 7(fc,a). It is defined via 

a ma 

where ft m (a) = p m /ptot- Note, however, that this parametrization assumes that the growth function is scale- 
independent, so it can only be used to probe for time-dependent modifications of growth. A significant deviation 
of the observed 7^ from its predicted value of 6/11 sa 0.55 would indicate a breakdown of ACDM. It is possible to use 
MGCAMB to constrain parameter 7^ by expressing p(k,a) in terms of 7^ as [30] 



^+2 + ^+7L^+7Un(a. 



(23) 

where the prime denotes differentiation with respect to In a, and the expression holds for a time-dependent 7^. 



III. THE PATCH 



CAMB [19, 20] is a publicly available code that evolves the Boltzmann and Einstein equations to calculate cosmo- 
logical observables such as the CMB and matter power spectra. The first version of MGCAMB [6] evolved the same 
equations as CAMB up to rcdshift z = 30, while at lower rcdshifts it used Eqs. (11) and (12), which neglect effects of 
relativistic species. Hence, as mentioned above, it was not suitable for describing modifications at early times {z > 30). 
This limitation is removed in the current version of MGCAMB, which uses Eqs. (13) and (14) valid for the entire 
redshift range and accounts for all the species. 

To implement the Newtonian gauge parametrization of (13) and (14) in CAMB, we first need to derive an equivalent 
system of equations in the synchronous gauge. We have [26] 

* = a + Ha , (24) 
$ = -Ha, (25) 

where a = (h + 6r])/2k 2 with h and 77 being the scalar metric potentials in the synchronous gauge defined in (3). Our 
Eqs. (13) and (14) then convert to 

k 2 (a + Ha) = -^p(k,a){ P A + 3(p + P)a}, (26) 
r 1 -Ha- 1 {a + Ha) = ^p{p + P)a, (27) 
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where k = 87rGa 2 . Following the notation in CAMB, we can define the perturbation to the expansion rate, Z, and 
shear, a* (not to be confused with the anisotropic stress a), as 

*-5- < 28) 

a* =ka. (29) 
In CAMB, the perturbed Einstein equations are used to evaluate these quantities: 3 

1 



k 1 ^ = kl-LZ - -up A , (30) 



c 2 

|fcV ~Z) = Kpqi (31) 



where pq = (p + P)v. As we are interested in testing modifications to GR, we can no longer use Einstein's equations, 
and instead need an alternative way of evaluating Z and a* based on Eqs. (13) and (14) combined with conservation 
and Boltzmann equations. We do this by deriving a and h = 2k 2 a — 6rj from known quantities, as explained below. 
To derive a, we substitute (26) into (27) and obtain 4 



a 



{r ? +g[ 7 pA + 3( 7 -l)(p + P)a]}/^, (32) 



where A is in synchronous gauge. Then, to get r), we first substitute a from (26) into (27) to write 



V = tia-!^T, (33) 



where 

r = 7A + 3(1 + w)a(j - 1) . (34) 

The time derivative of (33) gives 



r) = Ua + Ua- ^{2UT - 3H(1 + w)T + -T + f } . (35) 
2k z p 



Using the conservation equations in the above, we have the final expression for 77: 



i + |^(i + «0] 



+ A[p{ 1 -l)H- f n-jp} 
+ 3<r(l + w)(l - i)fi + 3<r(l + 10) [3wp(-y - 1)U - (7 - l)/i - p-y]} , (36) 



where V is 



V = k 2 + y7MP(l + w) ■ (37) 

With a, fj and h in hand, we can evaluate Z and <7* in Eqs. (28,29), which are then used in the remaining equations 
of CAMB. 



IV. A WORKED EXAMPLE: JOINT CONSTRAINTS ON MASSIVE NEUTRINOS AND MODIFIED 

GRAVITY 



To demonstrate the use of MGCAMB in conjunction with CosmoMC, we performed a joint analysis of CMB, ISW 
and SNe data to constrain a particular set of modified gravity (MG) parameters together with the neutrino mass. 
Massive neutrinos can modify the growth of structure on different scales and redshifts depending on their mass. This 



3 Note that the corresponding equations in "CAMB Notes" are written in a different notation: their h is equivalent to h/6 in our notation 

4 When evolving the equations of motion in CAMB, r] is known from the previous time step 
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modification can, in principle, be degenerate with the effects of MG on the overlapping scales and rcdshifts. To study 
the degeneracy between MG parameters and basic cosmological parameters, including the neutrino mass, we chose 
two specific parametrized MG models, namely, the f(R) and Yukawa-type models introduced in [7] (see Eq. (18)). 
Cosmological constraints on f(R) models were also studied in [42]. 

On linear scales, modified gravity models, e.g. f(R), generally predict enhanced growth within the Compton wave- 
length due to the fifth force, while massive neutrinos can damp the structure within its free-streaming scale, which 
is determined by the neutrino mass. Therefore, it is interesting to study the degeneracy between the MG parameters 
and the neutrino mass. The new version of MGCAMB released with this paper is just the tool needed to perform this 
analysis. Fig 1 shows the CMB TT power spectra for ACDM, an " f(R) model" described by Eq. (18) with B Q = 0.5 
and the case of massive neutrinos with /„ = Q u /n C DM = 0.05 corresponding to ^2,m v ~ 0.5eV. We attempt to 
constrain four models: the f(R) and Yukawa-type, described by (18), each with and without the presence of massive 
neutrinos. In all the cases, we vary the 7 primary parameters as explained in Table I. In addition, we vary /„ for the 
cases with neutrinos, Bq in the case of f(R) and {Bq, (3i, s} for the Yukawa model. 
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FIG. 1: CMB TT power spectrum for ACDM (black solid), f(R) model with B = 0.5 (red dashed) and ACDM+massive 
neutrino models with f v = 0.05 (green dash-dotted). 

We utilized MGCAMB to calculate the CMB angular spectrum and the CMB/matter cross-correlation spectrum, as 
well as the SNe and BAO, and used a modified version of CosmoMC to fit the models to data. We used the temperature- 
temperature (TT) and temperature-polarization (TE) power spectra from the WMAP seven year observation [43], 
and the ISW-large scale structure (LSS) cross-correlation data by Ho et al [25] (see Appendix B for the details of 
the implementation). We also use the SNe data (SDSS compilation) [44], HST data from [45] and BAO data from 
[46] to further constrain the background expansion history. In addition, we put a top hat prior of [10, 20] Giga years 
on the age of the universe. Given the observational data, we use CosmoMC to sample the parameter space using the 
Metropolis-Hastings algorithm. Table I shows the parameters we used for sampling and the assumed priors. We run 
several chains, and obtained the constraints after the chains converge perfectly. 

Figs. 2 and 3 show the 1-D posterior distributions, and the 2-D contour plots of the cosmological and MG parameters 
for the f(R) and Yukawa-type models with (left panels) and without (right panels) massive neutrinos. We did not 
see a significant correlation between the neutrino mass and the MG parameters. This is because the MG models we 
have considered primarily affect the CMB spectrum via the ISW effect, which is relevant on large scales, or small £. 
On the other hand, the effect of (small) neutrino mass on the CMB spectrum is quite subtle and the constraint comes 
primarily from £ around the acoustic peaks. 

The correlation between the neutrino mass and MG would be more prominent if we considered their effect on the 
matter power spectrum data P(k). However, adding the P(k) information is non-trivial for the following reason. The 
MCMC module of Ho et al [25] is a compilation of clustering data from several surveys, and P(k) from each data set 
is used in determining the bias, which is then used to determine the cross-correlation of clustering with CMB. So, in 
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Parameter 


Explanation 


range (min, max) 




Primary parameters 


f(R) 


Yukawa-type 


■a* 

T 

ln(10 10 ^) 

A sz 
n s 


physical baryon density; uib = h 2 Q.b 

physical CDM density; co c = h 2 fl c 

sound horizon angle; ■&„ — 100 ■ r s (z*)/ Da(z*) 

optical depth to reionisation 

As is the scalar primordial amplitude 

amplitude of the SZ template for WMAP and ACBAR 

spectral index of primordial perturbations; n, — 1 = d In P/d In k 


(0.005,0.100) 
(0.01,0.99) 
(0.5, 10.0) 
(0.01,0.80) 
(2.7,4.0) 

(0,2) 
(0.5,1.5) 




Neutrino parameters 


f(R) 


Yukawa-type 


u 


fraction of dark matter density as massive neutrinos 




(0.0,0.1) 




MG parameters 


f(R) 


Yukawa-type 


B 

Pi 
s 


present lengthscale of the theory (in units of the horizon scale) 
coupling 

time evolution of the scalaron mass 


(0,1) 
4/3 
4 


(0,1) 
(0.001,2) 

(1,4) 



TABLE I: List of the parameters used in the Monte Carlo sampling. The ranges of the fiat priors are given if a parameter is 
varied, or the value is given if the parameter was fixed. 

effect, P(k) is used, but only in determining the bias, and not in constraining cosmological parameters. Same method 
was used in [15] and in [47]. In principle, it should be possible to include P(k) into the data while properly accounting 
for the covariance with cross-correlation, but this task is outside the scope of this paper. 

For the f(R) model, we find that the Bq parameter is constrained to be Bq < 0.4 (95% C.L.), which corresponds to 
\\ < 1900 Mpc/h. The constraint on B is practically unchanged after marginalizing over the neutrino mass. For the 
Yukawa model, we find that 0.7 < (3i < 1.7 (95% C. L.), while Bo and s are very weakly constrained. In both cases, 
we find f v < 0.05 at 95% C.L., implying J2 m v ~ 0.5 Our constraints on the MG parameters are consistent with 
those presented in Refs. [7, 15]. The analysis of Rcf. [7] is based on a different compilation of ISW datasets - we used 
Ho et al. [25] , while they used the compilation of Giannantonio et al. [47] . Our results are also in good agreement with 
those of Ref. [15], where the PPF framework [3, 40] employing a different set of equations was used. When the Ho et 
al. [25] datasct was considered in [15], it lead to a constraint of Bq < 0.4 (95% C.L.) on f(R) models. We hope that 
including the weak lensing and peculiar velocity data [14] and the Cosmic Mach Number (CMN) [48] will improve 
the constraints, and it is straightforward to add additional observables given the tool we developed in this work. 




FIG. 2: Marginalized posterior distribution for the f(R) model parameters and 2-D contour plots showing the ranges of and 
correlations between parameters of interest and the 68 and 95% confidence limits, left: with massive neutrinos, right: without 
massive neutrinos. 
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FIG. 3: Same as in Fig 2 but for the Chameleon-type model considered in the paper. 



V. SUMMARY 

We have introduced a patch to CAMB and CosmoMC, called MGCAMB, that makes it possible to evaluate 
cosmological observables using a parametrized modification of linear Einstein equations. The motivation for this 
paramctrization in terms of functions fj,(a,k) and 7(0, k), its consistency, and the relation to other models in the 
literature has been discussed in [30]. Its main value is in allowing for mo del- independent way of looking for depar- 
tures from GR [8]. However, it is also possible to derive approximate expressions for fj,(a,k) and j(a, k) in terms of 
fundamental parameters of particular theories [7, 31]. 

We have demonstrated the advantages afforded by this new version of MGCAMB by obtaining joint constraints on 
the neutrino mass and parameters of two models of modified gravity previously considered in [7]. We find that GR 
remains a good fit in all cases. In the case of f(R) we find that the Bo parameter is constrained to be B$ < 0.4 (95% 
C.L.), which corresponds to \\ < 1900 Mpc/h. In the case of Yukawa we find that 0.7 < pi < 1.7 (95% C. L.), while 
Bq and s are very weakly constrained. In both models there is little degeneracy between the neutrino mass and the 
MG parameters, with the neutrino fraction constrained to be /„ < 0.05 at 95% C.L.. 

We have made the MGCAMB patch public at http://www.sfu.ca/~aha25/MGCAMB.html and will update it regu- 
larly to keep it compatible with the latest distributions of CAMB. 
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Appendix A: Alterations to the CMB source function in CAMB. 

In the notation of [49], the CMB temperature angular spectrum Cj T can be written as 

Cj T = (4tt) 2 / k 2 dkP(k)[Af (fc)] 2 , 



where P(k) is the primordial power spectrum and 

Aj(k) = 



d,TS T (k,T)ji(kT) , 



(Al) 



(A2) 



where r is the conformal time, Sx(k, r) is the source term, and je(x) are spherical Bcsscl functions. In GR, the source 
term in terms of the synchronous gauge variables is [49] 



S T (k,r) = g A 



n . v b n 3n 

to -i- 2a + — + — + — T 
k 4 4fc 2 



+ e K (ri + a)+g [ a+-+ n 



v b 3n 

I + 2k 2 



3gU 

4fc 2 



(A3) 



where k is the optical depth, g is the visibility function, II = + Af + Aq and Aj(Af) are the £'th moments of 
A T (A P ) in term of Legendre polynomials [49]. 

In the default version of CAMB, the perturbed Einstein equations are used in the evaluation of the source Sr(k, r). 
Since we have introduced modifications to the Poisson (26) and anisotropy (27) equations, it was important to avoid 
assuming GR in the expression for the source. The MGCAMB patch replaces the default CMB source in CAMB with 
the one given by Eq. (A3). In addition, the evaluation of the ISW term, e~ K (i) + a), in (A3) is modified. Namely, it 
was originally given by 



1] + a= 2kt 



while the MGCAMB patch changes it to 



2(pA + P A) + l{p + P)& + |(p + P)a 



( 7 + l)(pA + P A) + ^(p + P)& + 7 ^(p + P)a 



■7M 



(pA) + -{p + P)a 



(A4) 



(A5) 



where Eqs. (32-35) are used. It can be easily checked that for \x = 7 = 1, we get back the ISW term in GR. 



Appendix B: Alterations to the CMB-LSS cross-correlation patch of Ho et al. 

To use the ISW data of Ho et al. [25], we need to calculate the theoretical CMB-LSS cross-correlation angular 
spectrum, Cf. Below, we explain how it is done with MGCAMB. 
In [25], Cf T is written as 5 

Cf = - ( k 2 dkP(k)g e (k)T e (k) , (Bl) 



with 

9t(k) = I dzbi(z)IL(z)D(z)ji(k X (z)) , (B2) 



5 Note that the different pre-factor, (2/w) in Eq. (Bf) rather than (4-7r) 2 in (Eq. (Af)), is due to different Fourier conventions [50]. 
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and 



T e (k) = 3^0 m T CMB x J dz^j- [D(z)(l + z)] Je{k * {z)) , (B3) 

where II(z) is the normalized selection function, Tomb is temperature of CMB today, x( z ) is the comoving distance to 
rcdshift z, bi(z) is the bias factor, integration over r is replaced with integration over redshift and the growth factor 
D(z) is defined as 

S(k,z) D(z) 

6(k,0) D(0) ■ 1 ' 

Expression (Bl) can be simplified for £ > 10 by working in the flat-sky approximation and substituting k = (£ + 
1/2) /x(z), which gives [25] 

Cf = Zffm* X / **K*W*)^DM±[DizXl + z)]P (<±^) . (B5) 

In a general theory of gravity, the growth factor can have a complicated form depending both on scale and time. 
With Eqs. (26) and (27), it is straight-forward to show that at late times, when the anisotropic stress due to relativistic 
species is negligible, Eq. (B5) is modified to 

cf = Zffiw x / ^)n(^o(^P(,)(i + ,)MM)(i^(M))]p(^) , (as) 

where D(z) is still defined as Eq. (B4) and we use MGCAMB to calculate D(z) from Eq. (B4) at each redshift and 
also the modified matter power spectrum today, P((£ + l/2)/x). The derivative of the growth factor with respect to 
redshift can be calculated using Eqs. (8), (26) and the corresponding expressions for /J,(k,z) and 7(fc, z). 
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